Juvenile Pacific cod were tagged, acclimated to laboratory conditions, then gradually acclimated to four experimental temperatures (0, 5, 9, 16). After 6 weeks in treatments, fish were sacrificed, liver and fin tissues were collected for RNASeq and lcWGS, respectively, and measurements taken of fish length, wet weight, and liver weight.

On 11/21/22 160 tagged fish entered experimental tanks for acclimation, 40 per treatment (n=10 fish / tank, 4 tanks/treatment). On 12/28/22 temperatures began to slowly increase until all target experimental temperatures (0C, 5C, 9C, 16C) were reached on 1/8/23. Of the 40 fish per treatment, all survived the 0C, 5C, and 9C treatments. Four fish died in the 16degC treatment on 1/22/23 (tank 1), 1/24/23 (tank 2), 1/29/23 (tank 4), and 2/1/23 (tank 2). All survivors were euthanized/sampled on 2/8 & 2/9. Standard length and wet weight were 1) collected before tank acclimation, 2) before experimental temperature, and 3) at treatment termination so that growth and condition could be assessed prior to and during experimental treatments. Liver wet weight and tissues for sequencing were collected at treatment termination.

In this notebook I look at effects of temperature on distributions of growth rate (length, weight) body condition index (wet weight / length), and liver condition (i.e. hepato-somatic index, liver wet weight / whole body wet weight).

Load libraries and source scripts

list.of.packages <- c("tidyverse", "readxl", "janitor", "purrr", "ggpubr", "googlesheets4", "plotly", "lubridate")

# Load all required libraries
all.packages <- c(list.of.packages)
lapply(all.packages, FUN = function(X) {
  do.call("require", list(X))
})

`%!in%` = Negate(`%in%`)

Load data

load(file = "../rnaseq/aligned/counts.ts")

phenotypes <- read_excel("../data/Pcod Temp Growth experiment 2022-23 DATA.xlsx", sheet = "AllData") %>% clean_names() %>%
  mutate_at(c("tank", "temperature", "microchip_id", "dissection_date", "genetic_sampling_count"), factor) %>%
  dplyr::rename("sl_final"="sl_mm", "wwt_final"="whole_body_ww_g") %>% 
  mutate(growth.sl.accl=sl_12272022-sl_11212022,
         growth.sl.trt=round(sl_final-sl_12272022, digits = 2),
         growth.wwt.trt=wwt_final-wwt_12272022,
         growth.wwt.accl=round(wwt_12272022-wwt_11212022, digits = 2),
         ci=wwt_final/sl_final,
         hsi=total_liver_ww_mg/(wwt_final*1000),
         mort=as.factor(case_when(is.na(mort_date) ~ "No", TRUE ~ "Yes")),
         rna=case_when(genetic_sampling_count %in% gsub("sample_", "", rownames(counts.ts)) ~ "Yes", TRUE ~ "No"),
         measure_date=case_when(is.na(mort_date) ~ ymd(dissection_date), TRUE ~ mort_date)) %>%
  mutate(days_growth = as.numeric(difftime(ymd(measure_date), ymd("2022-12-27"), units = "days")))

phenotypes %>% head()
## # A tibble: 6 × 32
##   microchip_id sl_11212022 wwt_11212022 tank  temperature sl_12272022
##   <fct>              <dbl>        <dbl> <fct> <fct>             <dbl>
## 1 9397                  86         7.38 12    0                    93
## 2 9467                  92         7.98 12    0                    97
## 3 5059                 102        11.6  12    0                   113
## 4 9461                  82         6.15 12    0                    88
## 5 9560                  95         9.59 12    0                   101
## 6 7386                  95         9.23 12    0                   108
## # ℹ 26 more variables: wwt_12272022 <dbl>, mort_date <dttm>,
## #   dissection_date <fct>, sl_final <dbl>, wwt_final <dbl>,
## #   total_liver_ww_mg <dbl>, liverfor_lipids_ww_mg <dbl>,
## #   muscle_w_wfor_lipids_mg <dbl>, genetic_sampling_count <fct>,
## #   dissection_comments <chr>, molecular_notes <chr>, fin <chr>, liver <chr>,
## #   blood <chr>, spleen <chr>, gill <chr>, growth.sl.accl <dbl>,
## #   growth.sl.trt <dbl>, growth.wwt.trt <dbl>, growth.wwt.accl <dbl>, …
# Wet weight over time
phenotypes %>%
  group_by(temperature) %>%
  dplyr::summarise(sl_mean.1=mean(sl_11212022),sl_sd.1=sd(sl_11212022), 
            sl_mean.2=mean(sl_12272022),sl_sd.2=sd(sl_12272022),
            sl_mean.3=mean(sl_final),sl_sd.3=sd(sl_final),
            
            wwt_mean.1=mean(wwt_11212022),wwt_sd.1=sd(wwt_11212022), 
            wwt_mean.2=mean(wwt_12272022),wwt_sd.2=sd(wwt_12272022),
            wwt_mean.3=mean(wwt_final),wwt_sd.3=sd(wwt_final)) %>%
  pivot_longer(cols = -temperature) %>%
  separate(name, sep = "\\.", into = c("metric", "time")) %>%
  mutate_at(c("metric"), factor) %>% mutate(time=as.numeric(time)) %>%
  filter(metric=="wwt_mean") %>%
  ggplot() + geom_line(aes(x=time, y=value, color=temperature), cex=1.5) + theme_minimal() +
  scale_color_manual(name="Temperature", 
                       values=c("0"="#2c7bb6", "5"="#abd9e9", "9"="#fdae61", "16"="#d7191c")) + 
  ggtitle("Wet weight over time") + xlab(NULL) + ylab("Wet weight") +
  scale_x_continuous(breaks = c(1, 2, 3), labels = c("Pre-Acclimation", "Pre-Treatment", "Treatment\ntermination"))

### Standard length by temperature over time - mean and SD
phenotypes %>%
  group_by(temperature) %>%
  dplyr::summarise(sl_mean.1=mean(sl_11212022),sl_sd.1=sd(sl_11212022), 
            sl_mean.2=mean(sl_12272022),sl_sd.2=sd(sl_12272022),
            sl_mean.3=mean(sl_final),sl_sd.3=sd(sl_final),
            
            wwt_mean.1=mean(wwt_11212022),wwt_sd.1=sd(wwt_11212022), 
            wwt_mean.2=mean(wwt_12272022),wwt_sd.2=sd(wwt_12272022),
            wwt_mean.3=mean(wwt_final),wwt_sd.3=sd(wwt_final)) %>%
  pivot_longer(cols = -temperature) %>%
  separate(name, sep = "\\.", into = c("metric", "time")) %>%
  mutate_at(c("metric"), factor) %>% mutate(time=as.numeric(time)) %>%
  filter(metric=="sl_mean") %>%
  ggplot() + geom_line(aes(x=time, y=value, color=temperature), cex=1.5) + theme_minimal() +
  scale_color_manual(name="Temperature", 
                       values=c("0"="#2c7bb6", "5"="#abd9e9", "9"="#fdae61", "16"="#d7191c")) + 
  ggtitle("Standard length over time") + xlab(NULL) + ylab("Standard length") +
  scale_x_continuous(breaks = c(1, 2, 3), labels = c("Pre-Acclimation", "Pre-Treatment", "Treatment\ntermination"))

phenotypes %>% 
    ggplot(aes(y=growth.sl.accl/days_growth, x=temperature, fill=temperature, shape=mort, alpha=mort)) + 
  geom_violin(trim = F) + geom_point(position = position_jitterdodge(jitter.width = 1)) + theme_minimal() +
    ggtitle("Growth rate while acclimating, standard length (mm/day)") +
  scale_alpha_manual(values=c(0.75,0.5)) +
  scale_x_discrete(drop=T) + #Do drop empty factors for temperature contrasts
  ylab("mm / day") +
  scale_fill_manual(name="Temperature", values=c("0"="#2c7bb6", "5"="#abd9e9", "9"="#fdae61", "16"="#d7191c")) +
  xlab(NULL) + ylab(NULL)

# phenotypes %>% filter(mort=="No") %>%
#     ggplot(aes(y=growth.sl.accl/days_growth, x=temperature, fill=temperature)) + 
#   geom_violin() + geom_jitter(width=0.25) + theme_minimal() +
#     ggtitle("Growth rate while acclimating, standard length (mm/day)") +
#   scale_fill_manual(name="Temperature", 
#                        values=c("0"="#2c7bb6", "5"="#abd9e9", "9"="#fdae61", "16"="#d7191c")) + 
#   xlab(NULL) + ylab(NULL)

Growth while acclimating, standard length statistics

phenotypes %>% filter(mort=="No") %>% group_by(temperature) %>% dplyr::summarise(mean=mean(growth.sl.accl), sd=sd(growth.sl.accl))
## # A tibble: 4 × 3
##   temperature  mean    sd
##   <fct>       <dbl> <dbl>
## 1 0            8.7   2.82
## 2 5            8.43  3.34
## 3 9            8     2.42
## 4 16           9.78  2.75
summary(aov(growth.sl.accl ~ temperature, phenotypes %>% filter(mort=="No")))
##              Df Sum Sq Mean Sq F value Pr(>F)  
## temperature   3   64.4  21.481   2.637 0.0518 .
## Residuals   152 1238.4   8.147                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
TukeyHSD(aov(growth.sl.accl ~ temperature, phenotypes %>% filter(mort=="No")))
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = growth.sl.accl ~ temperature, data = phenotypes %>% filter(mort == "No"))
## 
## $temperature
##           diff         lwr       upr     p adj
## 5-0  -0.275000 -1.93297778 1.3829778 0.9730922
## 9-0  -0.700000 -2.35797778 0.9579778 0.6920681
## 16-0  1.077778 -0.62563246 2.7811880 0.3574449
## 9-5  -0.425000 -2.08297778 1.2329778 0.9097552
## 16-5  1.352778 -0.35063246 3.0561880 0.1700574
## 16-9  1.777778  0.07436754 3.4811880 0.0371523
phenotypes %>% 
    ggplot(aes(y=growth.sl.trt/days_growth, x=temperature, fill=temperature, shape=mort, alpha=mort)) + 
  geom_violin() + geom_point(position = position_jitterdodge(jitter.width = 1)) + theme_minimal() +
    ggtitle("Growth rate in treatments, standard length (mm/day)") +
  scale_alpha_manual(values=c(0.75,0.5)) +
  ylab("mm / day") +
  scale_fill_manual(name="Temperature", values=c("0"="#2c7bb6", "5"="#abd9e9", "9"="#fdae61", "16"="#d7191c")) +
  xlab(NULL) + ylab(NULL)

# phenotypes %>% filter(mort=="No") %>%
#     ggplot(aes(y=growth.sl.trt, x=temperature, fill=temperature)) + 
#   geom_violin() + geom_jitter(width=0.25) + theme_minimal() +
#     ggtitle("Growth in treatments, standard length") +
#   scale_fill_manual(name="Temperature", 
#                        values=c("0"="#2c7bb6", "5"="#abd9e9", "9"="#fdae61", "16"="#d7191c"))

Growth in treatments, standard length statistics

summary(aov(growth.sl.trt ~ temperature, phenotypes %>% filter(mort=="No")))
##              Df Sum Sq Mean Sq F value Pr(>F)    
## temperature   3   2456   818.8   87.06 <2e-16 ***
## Residuals   152   1430     9.4                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
TukeyHSD(aov(growth.sl.trt ~ temperature, phenotypes %>% filter(mort=="No")))
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = growth.sl.trt ~ temperature, data = phenotypes %>% filter(mort == "No"))
## 
## $temperature
##           diff       lwr        upr     p adj
## 5-0   4.875000  3.093614  6.6563863 0.0000000
## 9-0  10.275000  8.493614 12.0563863 0.0000000
## 16-0  8.569444  6.739244 10.3996449 0.0000000
## 9-5   5.400000  3.618614  7.1813863 0.0000000
## 16-5  3.694444  1.864244  5.5246449 0.0000031
## 16-9 -1.705556 -3.535756  0.1246449 0.0774644
phenotypes %>% filter(mort=="No") %>% group_by(temperature) %>% dplyr::summarise(mean=mean(growth.sl.trt), sd=sd(growth.sl.trt))
## # A tibble: 4 × 3
##   temperature  mean    sd
##   <fct>       <dbl> <dbl>
## 1 0            4.12  1.92
## 2 5            9     2.53
## 3 9           14.4   3.51
## 4 16          12.7   3.98
phenotypes %>% 
    ggplot(aes(y=growth.wwt.accl/days_growth, x=temperature, fill=temperature, shape=mort, alpha=mort)) + 
  geom_violin() + geom_point(position = position_jitterdodge(jitter.width = 1)) + theme_minimal() +
    ggtitle("Growth rate while acclimating, wet weight (g/day)") +
  scale_alpha_manual(values=c(0.75,0.5)) +
  ylab("g / day") +
  scale_fill_manual(name="Temperature", values=c("0"="#2c7bb6", "5"="#abd9e9", "9"="#fdae61", "16"="#d7191c")) +
  xlab(NULL) 

# phenotypes %>% filter(mort=="No") %>%
#     ggplot(aes(y=growth.wwt.accl, x=temperature, fill=temperature)) + 
#   geom_violin() + geom_jitter(width=0.25) + theme_minimal() +
#     ggtitle("Growth while acclimating, wet weight") +
#   scale_fill_manual(name="Temperature", 
#                        values=c("0"="#2c7bb6", "5"="#abd9e9", "9"="#fdae61", "16"="#d7191c"))

Growth while acclimating, wet weight - stats

summary(aov(growth.wwt.accl ~ temperature, phenotypes %>% filter(mort=="No")))
##              Df Sum Sq Mean Sq F value Pr(>F)
## temperature   3   3.64  1.2122    1.85   0.14
## Residuals   152  99.59  0.6552
TukeyHSD(aov(growth.wwt.accl ~ temperature, phenotypes %>% filter(mort=="No")))
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = growth.wwt.accl ~ temperature, data = phenotypes %>% filter(mort == "No"))
## 
## $temperature
##             diff        lwr       upr     p adj
## 5-0  -0.23050000 -0.7006701 0.2396701 0.5811808
## 9-0  -0.32350000 -0.7936701 0.1466701 0.2834053
## 16-0  0.04158333 -0.4414705 0.5246372 0.9960401
## 9-5  -0.09300000 -0.5631701 0.3771701 0.9556859
## 16-5  0.27208333 -0.2109705 0.7551372 0.4623147
## 16-9  0.36508333 -0.1179705 0.8481372 0.2065730
phenotypes %>% filter(mort=="No") %>% group_by(temperature) %>% dplyr::summarise(mean=mean(growth.wwt.accl), sd=sd(growth.wwt.accl))
## # A tibble: 4 × 3
##   temperature  mean    sd
##   <fct>       <dbl> <dbl>
## 1 0            2.45 0.994
## 2 5            2.22 0.682
## 3 9            2.13 0.811
## 4 16           2.50 0.703
phenotypes %>% 
    ggplot(aes(y=growth.wwt.trt/days_growth, x=temperature, fill=temperature, shape=mort, alpha=mort)) + 
  geom_violin() + geom_point(position = position_jitterdodge(jitter.width = 1)) + theme_minimal() +
    ggtitle("Growth rate in treatments, wet weight (g/day)") +
  scale_alpha_manual(values=c(0.75,0.5)) +
  ylab("g / day") +
  scale_fill_manual(name="Temperature", values=c("0"="#2c7bb6", "5"="#abd9e9", "9"="#fdae61", "16"="#d7191c"))

# phenotypes %>% filter(mort=="No") %>%
#     ggplot(aes(y=growth.wwt.trt, x=temperature, fill=temperature)) + 
#   geom_violin() + geom_jitter(width=0.25) + theme_minimal() +
#     ggtitle("Growth in treatments, wet weight") +
#   scale_fill_manual(name="Temperature", 
#                        values=c("0"="#2c7bb6", "5"="#abd9e9", "9"="#fdae61", "16"="#d7191c"))

Growth in treatments, wet weight - stats

summary(aov(growth.wwt.trt ~ temperature, phenotypes %>% filter(mort=="No")))
##              Df Sum Sq Mean Sq F value Pr(>F)    
## temperature   3  283.3   94.44   92.34 <2e-16 ***
## Residuals   152  155.4    1.02                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
TukeyHSD(aov(growth.wwt.trt ~ temperature, phenotypes %>% filter(mort=="No")))
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = growth.wwt.trt ~ temperature, data = phenotypes %>% filter(mort == "No"))
## 
## $temperature
##            diff        lwr       upr     p adj
## 5-0   1.4110000  0.8235964 1.9984036 0.0000000
## 9-0   3.3832500  2.7958464 3.9706536 0.0000000
## 16-0  2.9721389  2.3686391 3.5756387 0.0000000
## 9-5   1.9722500  1.3848464 2.5596536 0.0000000
## 16-5  1.5611389  0.9576391 2.1646387 0.0000000
## 16-9 -0.4111111 -1.0146109 0.1923887 0.2920273
phenotypes %>% filter(mort=="No") %>% group_by(temperature) %>% dplyr::summarise(mean=mean(growth.wwt.trt), sd=sd(growth.wwt.trt))
## # A tibble: 4 × 3
##   temperature   mean    sd
##   <fct>        <dbl> <dbl>
## 1 0           0.0518 0.526
## 2 5           1.46   0.770
## 3 9           3.44   1.25 
## 4 16          3.02   1.32
phenotypes %>% 
    ggplot(aes(y=ci, x=interaction(mort, temperature), fill=temperature, shape=mort, alpha=mort, label=genetic_sampling_count)) + 
  geom_violin() + geom_jitter(width=0.25) + theme_minimal() +
    ggtitle("Body condition index\n(wet weight / standard length") +
  scale_alpha_manual(values=c(0.75,0.5)) +
  scale_fill_manual(name="Temperature", 
       values=c("0"="#2c7bb6", "5"="#abd9e9", "9"="#fdae61", "16"="#d7191c")) +
  theme(axis.text.x = element_blank(), axis.ticks.x = element_blank()) +
  xlab(NULL) + ylab("Wet weight / standard length")

# phenotypes %>% filter(mort=="No") %>%
#     ggplot(aes(y=ci, x=temperature, fill=temperature)) + 
#   geom_violin() + geom_jitter(width=0.25) + theme_minimal() +
#     ggtitle("Appr. body condition index\n(wet weight / standard length") +
#   scale_fill_manual(name="Temperature", 
#                        values=c("0"="#2c7bb6", "5"="#abd9e9", "9"="#fdae61", "16"="#d7191c"))

Body condition index at termination (WWT/SL)

summary(aov(ci ~ temperature, phenotypes %>% filter(mort=="No")))
##              Df  Sum Sq  Mean Sq F value   Pr(>F)    
## temperature   3 0.00965 0.003216   10.39 2.91e-06 ***
## Residuals   152 0.04702 0.000309                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
TukeyHSD(aov(ci ~ temperature, phenotypes %>% filter(mort=="No")))
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = ci ~ temperature, data = phenotypes %>% filter(mort == "No"))
## 
## $temperature
##             diff          lwr        upr     p adj
## 5-0  0.002714042 -0.007502273 0.01293036 0.9007512
## 9-0  0.014829986  0.004613671 0.02504630 0.0013176
## 16-0 0.018819590  0.008323324 0.02931586 0.0000408
## 9-5  0.012115944  0.001899629 0.02233226 0.0129654
## 16-5 0.016105547  0.005609282 0.02660181 0.0005974
## 16-9 0.003989604 -0.006506662 0.01448587 0.7568554
phenotypes %>% filter(mort=="No") %>% group_by(temperature) %>% dplyr::summarise(mean=mean(ci), sd=sd(ci))
## # A tibble: 4 × 3
##   temperature  mean     sd
##   <fct>       <dbl>  <dbl>
## 1 0           0.104 0.0168
## 2 5           0.107 0.0140
## 3 9           0.119 0.0214
## 4 16          0.123 0.0174
phenotypes %>%
    ggplot(aes(y=hsi, x=interaction(mort, temperature), fill=temperature, shape=mort, alpha=mort)) + 
  geom_violin() + geom_jitter(width=0.25) + theme_minimal() +
    ggtitle("Hepato-somatic index\n(liver wet weight / total wet weight)") +
  scale_alpha_manual(values=c(0.75,0.5)) +
  scale_fill_manual(name="Temperature", 
       values=c("0"="#2c7bb6", "5"="#abd9e9", "9"="#fdae61", "16"="#d7191c")) +
  theme(axis.text.x = element_blank(), axis.ticks.x = element_blank()) +
  xlab(NULL) + ylab("Liver wet weight / Total wet weight")

# phenotypes %>% filter(mort=="No") %>%
#     ggplot(aes(y=hsi, x=temperature, fill=temperature)) + 
#   geom_violin() + geom_jitter(width=0.25) + theme_minimal() +
#     ggtitle("Appr. hepato-somatic index\n(liver wet weight / total wet weight)") +
#   scale_fill_manual(name="Temperature", 
#                        values=c("0"="#2c7bb6", "5"="#abd9e9", "9"="#fdae61", "16"="#d7191c"))

Hepato-somatic index - stats

summary(aov(hsi ~ temperature, phenotypes %>% filter(mort=="No")))
##              Df    Sum Sq   Mean Sq F value   Pr(>F)    
## temperature   3 4.945e-09 1.648e-09   31.45 7.14e-16 ***
## Residuals   152 7.966e-09 5.240e-11                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
TukeyHSD(aov(hsi ~ temperature, phenotypes %>% filter(mort=="No")))
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = hsi ~ temperature, data = phenotypes %>% filter(mort == "No"))
## 
## $temperature
##               diff           lwr           upr     p adj
## 5-0  -2.335143e-06 -6.540247e-06  1.869961e-06 0.4749567
## 9-0  -8.358065e-06 -1.256317e-05 -4.152961e-06 0.0000045
## 16-0 -1.473728e-05 -1.905762e-05 -1.041695e-05 0.0000000
## 9-5  -6.022922e-06 -1.022803e-05 -1.817818e-06 0.0015764
## 16-5 -1.240214e-05 -1.672247e-05 -8.081807e-06 0.0000000
## 16-9 -6.379219e-06 -1.069955e-05 -2.058885e-06 0.0010420
phenotypes %>% filter(mort=="No") %>% group_by(temperature) %>% dplyr::summarise(mean=mean(hsi), sd=sd(hsi))
## # A tibble: 4 × 3
##   temperature      mean         sd
##   <fct>           <dbl>      <dbl>
## 1 0           0.0000361 0.00000674
## 2 5           0.0000337 0.00000877
## 3 9           0.0000277 0.00000704
## 4 16          0.0000213 0.00000600

Explore possible index of performance in warming: Condition Index * Hepatosomatic Index

Within treatments, I’d like to use WGCNA to identify co-expressed genes that are linearly related to a performance indicator. Here, I explore using HSI*CI as the performance indicator. In 16C, HSI was lower than control (9C) but CI was higher than control, so fish with high CI & HSI could be “high performers”, and expression patterns could provide molecular indicators of performance.

# Condition Index * Hepato-somatic Index 
phenotypes %>% 
    ggplot(aes(y=ci*hsi, x=interaction(mort, temperature), fill=temperature, shape=mort, alpha=mort, label=genetic_sampling_count)) + 
  geom_violin() + geom_jitter(width=0.25) + theme_minimal() +
    ggtitle("Possible index of performance (?)\n(Condition index * Hepatosomatic Index)") +
  scale_alpha_manual(values=c(0.75,0.5)) +
  scale_fill_manual(name="Temperature", 
     values=c("0"="#2c7bb6", "5"="#abd9e9", "9"="#fdae61", "16"="#d7191c")) +
  theme(axis.text.x = element_blank(), axis.ticks.x = element_blank()) +
  xlab(NULL) + ylab("CI * HSI")

Index of performance (CI * HSI) - stats

summary(aov(hsi*ci ~ temperature, phenotypes %>% filter(mort=="No")))
##              Df    Sum Sq   Mean Sq F value   Pr(>F)    
## temperature   3 2.921e-11 9.738e-12   9.545 8.18e-06 ***
## Residuals   152 1.551e-10 1.020e-12                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
TukeyHSD(aov(hsi*ci ~ temperature, phenotypes %>% filter(mort=="No")))
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = hsi * ci ~ temperature, data = phenotypes %>% filter(mort == "No"))
## 
## $temperature
##               diff           lwr           upr     p adj
## 5-0  -1.648471e-07 -7.515427e-07  4.218485e-07 0.8849743
## 9-0  -4.366931e-07 -1.023389e-06  1.500025e-07 0.2184555
## 16-0 -1.156764e-06 -1.759537e-06 -5.539917e-07 0.0000099
## 9-5  -2.718460e-07 -8.585416e-07  3.148496e-07 0.6254920
## 16-5 -9.919170e-07 -1.594689e-06 -3.891446e-07 0.0001958
## 16-9 -7.200711e-07 -1.322843e-06 -1.172986e-07 0.0121118
phenotypes %>% filter(mort=="No") %>% group_by(temperature) %>% dplyr::summarise(mean=mean(hsi*ci), sd=sd(hsi*ci))
## # A tibble: 4 × 3
##   temperature       mean          sd
##   <fct>            <dbl>       <dbl>
## 1 0           0.00000378 0.00000101 
## 2 5           0.00000362 0.00000106 
## 3 9           0.00000335 0.00000110 
## 4 16          0.00000263 0.000000828
#ggplotly(
  phenotypes %>% filter(mort=="No") %>%
    ggplot(aes(y=hsi, x=ci, color=temperature)) + 
  geom_point() + theme_minimal() +
    ggtitle("HSI / CI") +
  scale_color_manual(name="Temperature", 
                       values=c("0"="#2c7bb6", "5"="#abd9e9", "9"="#fdae61", "16"="#d7191c"))+
  geom_smooth(method = "lm")#)

Do we have RNASeq data from liver for the “Top Performers”?

ggplotly(phenotypes %>% filter(mort=="No") %>%
    ggplot(aes(y=hsi*ci, x=interaction(rna, temperature), fill=temperature, label=genetic_sampling_count, shape=rna)) + 
  geom_violin() + geom_jitter(width=0.25) + theme_minimal() +
    ggtitle("Performance Index (CI*HSI) by temperature, by RNASeq data or none") +
  scale_fill_manual(name="Temperature", 
                       values=c("0"="#2c7bb6", "5"="#abd9e9", "9"="#fdae61", "16"="#d7191c")))
# ggplotly(phenotypes %>% filter(rna=="Yes") %>%
#     ggplot(aes(y=hsi*ci, x=temperature, fill=temperature, label=genetic_sampling_count)) + 
#   geom_violin() + geom_jitter(width=0.25) + theme_minimal() +
#     ggtitle("Appr. hepato-somatic index\n(liver wet weight / total wet weight)") +      
#     scale_alpha_manual(values=c(0.75,0.5)) + ggtitle("Performance index by temperature, RNASeq samples only") +
#   scale_fill_manual(name="Temperature", 
#                        values=c("0"="#2c7bb6", "5"="#abd9e9", "9"="#fdae61", "16"="#d7191c")))
# 
# ggplotly(phenotypes %>% filter(rna=="Yes") %>%
#            left_join(haplos.allzp3, by = c("genetic_sampling_count"="sample")) %>%
#     ggplot(aes(y=hsi*ci, x=temperature, fill=temperature, label=genetic_sampling_count, shape=haplo.exons.lcwgs)) + 
#   geom_violin() + geom_jitter(width=0.25) + theme_minimal() +
#     ggtitle("Appr. hepato-somatic index\n(liver wet weight / total wet weight)") +
#   scale_fill_manual(name="Temperature", 
#                        values=c("0"="#2c7bb6", "5"="#abd9e9", "9"="#fdae61", "16"="#d7191c")))

Lipid data!

lipids <- read_excel("../data/Lipid class liver data_091324.xlsx") %>% clean_names() %>%
  mutate(sample_id=as.factor(sample_id))
## New names:
## • `TAG` -> `TAG...20`
## • `FFA` -> `FFA...21`
## • `ST` -> `ST...22`
## • `Polar Lipids` -> `Polar Lipids...23`
## • `Total` -> `Total...24`
## • `TAG` -> `TAG...28`
## • `FFA` -> `FFA...29`
## • `ST` -> `ST...30`
## • `Polar Lipids` -> `Polar Lipids...31`
## • `Total` -> `Total...32`
phenotypes2 <- phenotypes %>% 
  left_join(lipids %>% dplyr::select(sample_id, genetics_sample_number, tag_20:total_32), 
            by=c("microchip_id"="sample_id")) %>%
  mutate(pi=1000*growth.sl.trt*ci*hsi*total_32)
phenotypes2 %>% 
  filter(!is.na(total_32)) %>%
    ggplot(aes(y=tag_28, x=interaction(mort, temperature), fill=temperature, shape=mort, alpha=mort, label=genetic_sampling_count)) + 
  geom_violin() + geom_jitter(width=0.25) + theme_minimal() +
    ggtitle("Triglyceride content (ug/mg)") +
  scale_alpha_manual(values=c(0.75,0.5)) +
  scale_fill_manual(name="Temperature", 
     values=c("0"="#2c7bb6", "5"="#abd9e9", "9"="#fdae61", "16"="#d7191c")) +
  theme(axis.text.x = element_blank(), axis.ticks.x = element_blank()) +
  xlab(NULL) + ylab("TAG Content")
## Warning: Groups with fewer than two data points have been dropped.
## Warning: The following aesthetics were dropped during statistical transformation: label
## ℹ This can happen when ggplot fails to infer the correct grouping structure in
##   the data.
## ℹ Did you forget to specify a `group` aesthetic or to convert a numerical
##   variable into a factor?

phenotypes2 %>% 
  filter(!is.na(total_32)) %>%
    ggplot(aes(y=ffa_29, x=interaction(mort, temperature), fill=temperature, shape=mort, alpha=mort, label=genetic_sampling_count)) + 
  geom_violin() + geom_jitter(width=0.25) + theme_minimal() +
    ggtitle("Free fatty acid content (ug/mg)") +
  scale_alpha_manual(values=c(0.75,0.5)) +
  scale_fill_manual(name="Temperature", 
     values=c("0"="#2c7bb6", "5"="#abd9e9", "9"="#fdae61", "16"="#d7191c")) +
  theme(axis.text.x = element_blank(), axis.ticks.x = element_blank()) +
  xlab(NULL) + ylab("FFA Content")
## Warning: Groups with fewer than two data points have been dropped.
## The following aesthetics were dropped during statistical transformation: label
## ℹ This can happen when ggplot fails to infer the correct grouping structure in
##   the data.
## ℹ Did you forget to specify a `group` aesthetic or to convert a numerical
##   variable into a factor?

phenotypes2 %>% 
  filter(!is.na(total_32)) %>%
    ggplot(aes(y=st_30, x=interaction(mort, temperature), fill=temperature, shape=mort, alpha=mort, label=genetic_sampling_count)) + 
  geom_violin() + geom_jitter(width=0.25) + theme_minimal() +
    ggtitle("Sterol content (ug/mg)") +
  scale_alpha_manual(values=c(0.75,0.5)) +
  scale_fill_manual(name="Temperature", 
     values=c("0"="#2c7bb6", "5"="#abd9e9", "9"="#fdae61", "16"="#d7191c")) +
  theme(axis.text.x = element_blank(), axis.ticks.x = element_blank()) +
  xlab(NULL) + ylab("sterol Content")
## Warning: Groups with fewer than two data points have been dropped.
## The following aesthetics were dropped during statistical transformation: label
## ℹ This can happen when ggplot fails to infer the correct grouping structure in
##   the data.
## ℹ Did you forget to specify a `group` aesthetic or to convert a numerical
##   variable into a factor?

phenotypes2 %>% 
  filter(!is.na(total_32)) %>%
    ggplot(aes(y=polar_lipids_31, x=interaction(mort, temperature), fill=temperature, shape=mort, alpha=mort, label=genetic_sampling_count)) + 
  geom_violin() + geom_jitter(width=0.25) + theme_minimal() +
    ggtitle("Polar lipid content (ug/mg)") +
  scale_alpha_manual(values=c(0.75,0.5)) +
  scale_fill_manual(name="Temperature", 
     values=c("0"="#2c7bb6", "5"="#abd9e9", "9"="#fdae61", "16"="#d7191c")) +
  theme(axis.text.x = element_blank(), axis.ticks.x = element_blank()) +
  xlab(NULL) + ylab("Polar lipid Content")
## Warning: Groups with fewer than two data points have been dropped.
## The following aesthetics were dropped during statistical transformation: label
## ℹ This can happen when ggplot fails to infer the correct grouping structure in
##   the data.
## ℹ Did you forget to specify a `group` aesthetic or to convert a numerical
##   variable into a factor?

phenotypes2 %>% 
  filter(!is.na(total_32)) %>%
    ggplot(aes(y=total_32, x=interaction(mort, temperature), fill=temperature, shape=mort, alpha=mort, label=genetic_sampling_count)) + 
  geom_violin() + geom_jitter(width=0.25) + theme_minimal() +
    ggtitle("Total lipid content (ug/mg)") +
  scale_alpha_manual(values=c(0.75,0.5)) +
  scale_fill_manual(name="Temperature", 
     values=c("0"="#2c7bb6", "5"="#abd9e9", "9"="#fdae61", "16"="#d7191c")) +
  theme(axis.text.x = element_blank(), axis.ticks.x = element_blank()) +
  xlab(NULL) + ylab("Total lipid Content")
## Warning: Groups with fewer than two data points have been dropped.
## The following aesthetics were dropped during statistical transformation: label
## ℹ This can happen when ggplot fails to infer the correct grouping structure in
##   the data.
## ℹ Did you forget to specify a `group` aesthetic or to convert a numerical
##   variable into a factor?

# Condition Index * Hepato-somatic Index * Total Lipid
#ggplotly(

phenotypes2 %>%
  # First I scale each metric so they are between about 1-20
    mutate(hsi=hsi*1e5,
           ci=ci*100,
           total_32=total_32/100) %>%
  #then I calculate a performance index from all metrics (x100)
    mutate(pi=(growth.sl.trt*hsi*ci*total_32)/100) %>%

      ggplot(aes(y=pi, x=interaction(mort, temperature), fill=temperature, shape=mort, alpha=mort, label=genetic_sampling_count)) + 
  geom_violin() + geom_jitter(width=0.25) + theme_minimal() +
    ggtitle("Possible composite index of performance\n(Growth * Condition index * Hepatosomatic Index * Total Lipid)") +
  scale_alpha_manual(values=c(0.75,0.5)) +
  scale_fill_manual(name="Temperature", 
     values=c("0"="#2c7bb6", "5"="#abd9e9", "9"="#fdae61", "16"="#d7191c")) +
  theme(axis.text.x = element_blank(), axis.ticks.x = element_blank()) +
  xlab(NULL) + ylab("Growth + CI * HSI * Total Lipid * 1,000")

#)